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Abstract 

The Porto Oscillation Code (POSC) has been developed 
in 1995 and improved over the years, with the main goal of 
calculating linear adiabatic oscillations for models of solar- 
type stars. It has also been used to estimate the frequencies 
and eigenfunctions of stars from the pre-main sequence up 
to the sub-giant phase, having a mass between 0.8 and 4 
solar masses. 

The code solves the linearised perturbation equations of 
adiabatic pulsations for an equilibrium model using a sec- 
ond order numerical integration method. The possibility of 
using Richardson extrapolation is implemented. Several op- 
tions for the surface boundary condition can be used. In this 
work we briefly review the key ingredients of the calcula- 
tions, namely the equations, the numerical scheme and the 
output. 

Keywords stars: interiors stars; oscillations methods: nu- 
merical 



1 Introduction 

The Porto Oscillation Code (POSC) was initially developed 
in 1995 to obtain the frequencies of solar models and en- 
velopes. The first description of the code has been given in 
Monteiroi ri996). 

The objective of this paper is to present a summary on 
how POSC calculates the frequencies of oscillations for stel- 
lar models. The paper starts with the basic linear equations 
describing the oscillations and how these are formulated to 
be solved numerically. The boundary conditions used and 
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their implementation are also discussed as well as the ac- 
curacy of the calculations. We end by listing some of the 
output values provided by the code and some of the applica- 
tions where the results of the code have been used. 



2 Basic equations for linear perturbations 

Our objective here is to review the necessary equations for 
non-radial adiabatic oscillations of spherically sy mmetric 
non-ro tating stars. By following the work by lUnno et al 



(|l989') it is possible to start from the hydrodynamic equa- 
tions (continuity, Poisson and conservation of momentum 
equations), in order to obtain a set of equations describ- 
ing the radial dependence of the amplitude functions for 
small perturbations. These perturbations correspond to; P 
for pressure, <I> for gravitational potential, while ^ is the dis- 
placement. The solutions are writen as 



P(f,r,0,0)=PoW+P(r) Yr{e,<P)e"" , 
4>(f,r,0,0)=4>o(r) + 4>(r) Yr{d,<^)e"'\ 

drr Ur) drr 
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Where the equilibrium configuration of the stars is described 
by the functions; po (for density), Pq and <I>o. Here t is time, 
CO the frequency for the oscilating solutions, (0,0) the hor- 
izontal variables while r is radial distance and Yi"^{9,(p) the 
spherical harmonics characterized by the integer numbers I 
(mode degree) and m (azimutal order with m^—l, ..,0, ..,/). 

By considering an equation for adiabatic perturbations 
and after eliminating the dependence on the horizontal coor- 
dinates and time, the equations describing the radial ampli- 
tude of the small perturbations are obtained in the following 
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form; 
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These form the set of equations we need to solve to ob- 
tain the radial behaviour of linear adiabatic oscillations of 
spherically symmetric stars. 
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The equilibrium structure in these equations is also char- 
acterized by quantities as gravity go and sound speed cq; 



(2) 3 The equilibrium model 



In order to describe the reference/equilibrium model we 
consider the following dimensionles s functions of the equi- 
libriu m structure (as defined by IChristensen-Dalsgaard 
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where the derivate has been calculated at fixed entropy S. 
There are also two characteristic frequencies; the Lamb fre- 
quency Si and the buoyancy frequency No (also known as the 
Brunt-Vaissala frequency), corresponding to 
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If we consider the following dimensionless variables 
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the equations can be written as 
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where M and are respectively the total mass and radius of 
the star, while m, is the mass within a sphere of radius r. 
These 5 functions are the result of an evolution code, being 
necessary as the input of the oscillation code. 

The four first order differential equations for small am- 
plitudes can now be written simply as 
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dx 
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where we have introduced the reduced frequency 
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We may write these equations in a vectorial form by 
defining the matrix, with L^=l{l+l), 



(10) 
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The system of differential equations is then simply written 

as 



ax 

where the vector y has the components {y\^y2^y^^y^)■ 



(11) 



4 Boundary conditions 

To complete the required equations it is also necessary to 
define four boundary conditions. The solution is to be found 
by integrating the equations between the centre of the star 
(r=0 or x=0) and the top of the atmosphere {r>R or x>X). 
So, in fact we shall be establishing two boundary conditions 
at x=0 and other two at the surface. The result is an eigen- 
value problem with solutions existing for discrete values of 
a. These are the eigenvalues associated to the corresponding 
eigenfunctions, that must satisfy the boundary conditions. 

4. 1 At the centre 

Since we have four dependent variables, the interior bound- 
ary conditions correspond to fix the values for two of the 
dependent variables. The other two are then related to these. 

The boundary conditions have to guarantee that the so- 
lutions are regular in the singular point, x=0, of the dif- 
ferential equations. So, we start by determining the limit- 
ing behaviour of the fl,'s when x^O. Considering that for 
x=r//?^l we can write (the subscript "c" stands for the 
value atx=0); 



p ^ Pc ■, nir 
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These expressions determine the behaviour of the fl,'s 
near the centre as follows from the definitions (|7]i; 

An 

ai ^ — Pc , 02 ^ , r^Tic , 04 , as 3 , (13) 

where Pc^R^Pc/M. 

If we now replace these in the definition ( fTOl l of the ma- 
trix Mf, the problem is reduced to a simple system of differ- 
ential equations with constant coefficients. This is. 
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This has a general solution (non-zero) given by 
; 7=1,2,3,4, 
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and 



%A =1 %3. 
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These are the two boundary conditions at the centre: from 
the values of ™d it is possible to determine and 

4.2 At the atmosphere 

Two more boundary conditions need to be imposed at the 
top of the atmosphere. In a similar fashion to what has been 
done for the centre, we now need to established what is the 
limiting behaviour for the a's (the subscript "5" represents in 
the following the value at the top of the atmosphere, located 
at f'fi/R>l) forx— >Xs. 

There are different options for imposing a boundary con- 
dition at the top of the model (surface). The most commonly 
used one is to assume an isothermal atmosphere for which 
we have that. 



0\ ^ 1 . 02 ^ 02s 7 ^3 ~ Tlv 7 O4 ~ 02s , ^5 ~ . 



(18) 



In such an isothermal atmosphere the density decreases ex- 
ponentially with radius. This behaviour allows one to ap- 
proximate the actual value of 055 by zero. 
The set of equations is now written as. 
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where the matrix has been approximated using the ap- 
proximate values of a, (seefTSli: 
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By redoing the analysis presented in the previous subsec- 
tion, and using 



yy=x-' £^,,;x2' ;;=1,2,3,4, 

i=0 

it follows that 

{L^/a^ ~ 02s) %a2 + fl2s^s,03 



(21) 
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^5,04 = -(' + 1)^5.03 . 



(22) 
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When using the first expression one must be careful since 
the actual value for 7'/^ can be imaginary. If it happens the 
solution will have a propagating component at the boundary, 
which implies that the wave will be loosing energy at this 
boundary. This does not correspond to the type of solutions 
we are looking for (standing waves). Therefore we only 
consider eigenvalues that are real, corresponding to standing 
waves, i.e. solutions that are evanescent at the boundaries. 

Note that this imposes restrictions on the values the fre- 
quency (7 can have for possible modes of oscillation. Solu- 
tions are only calculated for 7 > 0. 

Other options for surface boundary conditions are pos- 
sible (and have been implemented in POSC). The simplest 
option is to impose full reflection at the top of the model. 
Such a condition is achieved by setting 8P—Q at x—x^, giv- 
ing that 

^^5.01=^^5.02 + ^5,03, (23) 

instead of the first expression in Eqs (l22T i. 
5 Calculation of the solutions 

In order to calculated the eigenvalues (frequencies of oscilla- 
tion) of a solar model POSC uses a simple numerical scheme 
to solve the set of equations (fTTT i with the boundary condi- 
tions (Tl\ and ( |22] | (or one of the other alternatives). In this 
Section we describe briefly how this is done. 

The actual expressions implemented in the code are ex- 
tracted from the basic dimensionless system of 4 differential 
equations; 



dy 
dx 



■y ■ 



(24) 



where the matrix £/ is given in Eq. ( fTOl l. 

The functions a; are as listed in Eq. (|7| and known on a 
mesh from the equilibrium model obtained from solving the 
stellar structure equations. We use the dimensionless fre- 
quency (7, related to the actual frequency of oscillation {(o). 
As discussed above, under the selected boundary conditions, 
solutions exist only for discrete values of O = O/,,. The 
mode order n is associated with the radial structure of the 
different eigenfunctions that exist for the same mode degree 
/ (the code considers spherical stars, and so the solutions are 
independent of the azimutal order m). 

These values and the corresponding solution y are what 
we are trying to find. The method we use consists in, given 
a value of the degree I, to determine the values(s) of a that 
give a continuous solution at some meeting point (defined 
below as Xf). This point is where we stop the integration 
up from the centre, and the integration down from the atmo- 
sphere. In other words we find the value(s) of a that have a 
global solution satisfying all our four boundary conditions. 



So what we do in fact is to iterate in a in order to find the 
values that give the zeros of a function measuring the fitting 
of outer and inner solutions at x^. 

5.1 Numerical variables 

Due to numerical control of errors and precision of the cal- 
culation we redefine the variables for different regions of the 
model. We do so by estimating which regions of the star are 
evanescent for a given frequency. The two points used here 
to define these regions are Xin and Xout- These depend on the 
model and the values of (£> and degree /, corresponding to 
the roots of the following equation. 



(25) 



The acoustic cutoff frequency (Oc used here is determined by 
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(26) 



For the inner (near the centre) evanescent region we re- 
define the variables according to 



2-1 



y- 



(27) 



Here, Xin is the transition point separating this inner re- 
gion from the zone where the default variables, as given in 
Eq. (l24l i. are used. 

For the outer evanescent region (surface layers), above 
x=Xout, we use instead 



Jout = ( — ) y- 

\Xoul J 



(28) 



The equations are integrated from x~Q to x=Xin deter- 
mining Jin. From there to a fitting point Xf (well within 
the oscillatory region) we calculate the solution using the 
equations for y. Note that the transition from one region to 
the other is quite natural considering our definitions and 
Jout of y. On the other hand we integrate inward from x=Xs 
to x=Xout using instead the equations for jout- From there, 
down to X/ we take again the equations for y. 

Resulting from these two integrations we have the two 
sets of values at x^Xf which are then continuous (after nor- 
malization). Since the system of equations is linear, this is 
so if and only if the value of a is an eigenvalue. At this point 
what we actually do is to iterate on o to find the zeros of the 
fitting determinant at xy^. 

5.2 Method of integration 

The method implemented to solve numerically the equations 
considered above is a shooting method using a second-order 
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differences representation of the equations. It consists in 
writing the differential equations relating the values at two 
mesh points, x„ and x„+i, as 



y(«+i)=y(n) + - 



(«) + g(«+i) 



, (29) 



where h„ =x„^i —x„ and with y{n ) =y{x„ ) . In order to replace 
the derivatives we use the different sets of differential equa- 
tions discussed above for the regions O^Xjn ^^xj^^Xmit ■ 

We also have to implement the boundary conditions. It 
is done by setting the values of yi and y^ at the boundaries 
(centre and surface) and to calculate the values of y2 and 
y4 (at both boundaries) from the relations constructed in the 
previous Section. Both linearly independent solutions are 
found by setting the central/atmospheric values of yi equal 
to one and y^ alternatively to one and to zero. The actual 
solution is a linear combination of these two (for the interior 
solution - up toxf, as well as for the external solution - down 
toxf). 

Since we are using a shooting method, from the values of 
these two solutions ("in" and "out") at Xf, we construct the 
matching matrix whose determinant has to be zero if <7 is an 
eigenvalue. So the task of finding an eigenvalue is reduced to 
finding the zero of the determinant for the fitting conditions 
atXf. 

We maximize the efficiency of the search for the eigen- 
values (zeros of the determinant) by using the fact that these 
values are separated approximately by 



1/2 



A(7t 



fGM\ 



2na 



(p-modes). 



co'dr 



(30) 



(g-modes). 



dr 



5.3 Accuracy of the results 

The actual accuracy of the final values of the frequencies are 
determined by several aspects of the calculation. As it would 
be expected, the accuracy of the results (frequencies) de- 
pends on the accuracy of the equilibrium model being used. 
Here we do no t address t his issue, refer r ing th e reader to 
Monteiro et"al] (12006.) and lLebreton et al.l (l2008b . 

But another aspect associated with the equilibrium model, 
and determining the precision of the calculated eigen- 
values, is the mesh on which the equilibrium model is 
given. To minimise this effect, before calculating the 
frequencies we produce a re-meshing of the equilibrium 
model. The actual details of the new mesh depends on 
the type of model and oscillation modes being calcu- 
lated. We use a receipt similar to the one discussed by 
Christensen-Dalsgaard and Berthomieul (Il99lh . This allows 



us to minimize the errors caused by having too few points 
where the eigenfunctions are expected to vary more strongly. 

Other aspect determining the accuracy of the eigenvalues 
is of course the numerical method used to integrate the four 
differential equations discussed above. In our case we have 
a second order scheme for the integration of the system of 
differential equations. To this we have also added the use of 
reduced dependent variables in the regions where the ampli- 
tudes of the eigenfunctions would be otherwise very small. 

Further to this the code also uses an extrapolation to im- 
prove the accuracy of the determination of each frequency. 
It is known as Richardson Extrapolation. This uses the fact 
that our second order integration has an error which varies 
with the inverse of the squared number of mesh points. Us- 
ing such a fact it can be written that the actual value of the 
eigenvalue is 



(t2 = 



a 



On - ■ 



1 



a-1 



with a = I — 



(31) 



where Gn is the result found for a mesh of points and ff^/ 
for a mesh of A^' points. The code uses, by default, N'^N/2, 
giving a~4. This extrapolation requires extra work but im- 
proves si gnificantly the a ccuracy of the numerical frequen- 
cies (see iMova et al ] E008.) . 
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v (/liHz) 

Fig. 1 Comparison of the frequencies obtained using an equilib- 
rium model provided with a different number of mesh points. The 
reference is for a 10k mesh. Only differences for frequencies with 
/ = 0, 1,2,3 and 100 /xHz < v < 3500 /iHz are shown. 

The behaviour of the frequencies with increasing number 
of mesh points is an internal checkpoint that allows one to 
identify where the actual frequencies are no longer affected 
by the precision of the integration scheme. Such a compar- 
ison is shown in Fig. [T] where the frequencies obtained for 
a solar model in a mesh of 10k points is compared with the 
frequencies obtained using the same model with 4k and 8k 
mesh points. We have also perfor med a detailed compar- 
ison with the results from ADIPLS (iChristensen-Dalsgaard 
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to have an external check on the computation. Further 
compariso ns of POSC with oth er codes has been performed 
recently bv lMova et al] hoOSi) . 

In general (when using a model in a mesh of ^6k points) 
POSC frequencies of oscillation for solar-type stars have an 
estimated numerical uncertainty below 0.001 /xHz. This is 
below the current observational errors of solar frequencies. 
Similar values are obtained for oscillation modes of stars of 
different masses and ages if the mesh is adequately adapted 
(in terms of number of points and its distribution) to the 
eigenfunctions being calculated [g or p modes). 



6 Output 

The main output of the code are the values of the frequencies 
of linear adiabatic oscillations of an equilibrium model of a 
star The calculation to be done is defined by an interval in 
frequency (cOa < CO < COh) and in mode degree (la < I < lb)- 

6.1 Mode classification 

The mode order of each eigenvalue is obtained using 
a method s imilar to the phase diagram as described by 
Unno et al.l (Il989h . It consists in counting the number of 
times the solution crosses the line yi =0 in the plane {yi ,y2)- 
If the cross is clockwise it counts as (-1) otherwise as (+1). 
When /=0 (radial modes) an additional (+1) cross is consid- 
ered. 

To a total negative counting of the crosses corresponds a 
g-mode while p-modes have positive counting results, with 
the number corresponding to the mode order The solutions 
corresponding to /-mode eigenvalues have a total of zero 
counts. 



7 Conclusion 

This work provides a brief description of POSC - the Porto 
Oscillation Code. This code has been developed mainly for 
calculating linear adiabatic oscillations of stellar models for 
stars similar to the Sun (in mass). The code is written in 
Fortran 7 7 and is modular. It is prepared to accept in- 
put models in the AMDL0 format. Tools are also available 
to convert almost any available stellar model output to the 
required format to be used by POSC. 

T he code has been applied to several cases, namely the 



Sun ([Monteiroll99"6 : Monteiro et a l. 1996; Monteiro and Thomp son^ 
2005^ and other stars ( Cunha et al.^2003i : iFer nandes and M onteiroJ 
2003) , including pre-main sequence models (Ruoppoet aO 
2OO7I) . It has also been used to produce the frequencies of 



referenc e grid s of stellar evolution models for asteroseis- 
mology (IMarques et al.i.2008.) . 
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work was supported in part by the European Helio- and 
Asteroseismology Network (HELAS), a major international 
collaboration funded by the European Commission (FP6), 
as well as by FCT and POCI2010 (FEDER) through projects 
POCI/CTE-AST/57610/2004 and POCI/V.5/B0094/2005. 



6.2 Mode inertia and eigenfunctions 

The eigenfunctions are provided in different formats (sev- 
eral normalisations and/or combinations) depending on 
what is required. These are obtained from y and correspond 
to combinations of the functions. 



^rir) ^,{r) P{r) <I>(r) 



^r{Ry URY p{Ry hr)' 



(32) 



The equilibrium structure is also used to calculate different 
normalizations of the eigenfunctions. 

The code provides, in addition to the mode parameters 
and frequencies, the mode inertia as given by. 



Eln 



An 



(33) 



'See the description of some file formats for stellar evolution models at 
|http : / /www ■ astro ■ up ■ pt /c orot /ntoolsTl 
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